ANALYSIS, SYNTHESIS, AND ESTIMATION OF 
FRACTAL-RATE STOCHASTIC POINT 

PROCESSES 



On 
On 



STEFAN THURNER, STEVEN B. LOWEN, MARKUS C. FEURSTEIN 

and CONOR HENEGHAN 

Department of Electrical & Computer Engineering, Boston University, Boston, MA 02215 



ft. HANS G. FEICHTINGER 

on 

o 

m 



Institut fur Mathematik, Universitdt Wien, A-1090 Vienna, Austria 

MALVIN C. TEICH Q 

; Departments of Electrical & Computer Engineering, Biomedical Engineering, and Physics, 

§ ■ Boston University, Boston, MA 02215 

O 

r-> 

Abstract 

o 

& : 

els for describing a broad range of diverse phenomena, including electron transport in amor- 



Fractal and fractal-rate stochastic point processes (FSPPs and FRSPPs) provide useful mod- 



phous semiconductors, computer-network traffic, and sequences of neuronal action potentials. 
5_i ■ A particularly useful statistic of these processes is the fractal exponent a, which may be esti- 



mated for any FSPP or FRSPP by using a variety of statistical methods. Simulated FSPPs 
and FRSPPs consistently exhibit bias in this fractal exponent, however, rendering the study 
and analysis of these processes non-trivial. In this paper, we examine the synthesis and 
estimation of FRSPPs by carrying out a systematic series of simulations for several different 
types of FRSPP over a range of design values for a. The discrepancy between the desired 
and achieved values of a is shown to arise from finite data size and from the character of the 
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point-process generation mechanism. In the context of point-process simulation, reduction 
of this discrepancy requires generating data sets with either a large number of points, or 
with low jitter in the generation of the points. In the context of fractal data analysis, the 
results presented here suggest caution when interpreting fractal exponents estimated from 
experimental data sets. 
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1 INTRODUCTION 



Some random phenomena occur at discrete times or locations, with the individual events 
largely identical, such as a sequence of neural action potentials. A stochastic point process |L| 
is a mathematical construction which represents these events as random points in a space. 
Fractal stochastic point processes exhibit scaling in all of the statistics considered in this 
paper; fractal-rate stochastic point processes do so in some of them. In this work we consider 
the simulation and estimation of fractal and fractal-rate stochastic point processes on a line, 
which model a variety of observed phenomena in the physical and biological sciences. This 
work provides an extension and generalization of an earlier paper along these lines ||. 

1.1 Mathematical Descriptions of Stochastic Point Processes 

Figure [I] shows several representations that are useful in the analysis of point processes. 
Figure |l](a) demonstrates a sample function of a point process as a series of impulses occurring 
at specified times t n . Since these impulses have vanishing width, they are most rigorously 
defined as the derivative of a well-defined counting process N(t) [Fig. |l|(b)], a monotonically 
increasing function of t, that augments by unity when an event occurs. Accordingly, the 
point process itself is properly written as dN(t), since it is only strictly defined within the 
context of an integral. 

The point process is completely described by the set of event times {t n }, or equivalently 
by the set of interevent intervals. However, the sequence of counts depicted in Fig. 0(c) also 
contains much information about the process. Here the time axis is divided into equally 
spaced contiguous counting windows of duration T sec to produce a sequence of counts 
{Zk}, where Zk = N[(k + 1)T] — N[kT] denotes the number of events in the kth window. As 
illustrated in Fig. |](d), this sequence forms a discrete-time random process of nonnegative 
integers. In general, information is lost in forming the sequence of counts, although for a 
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regular point process the amount lost can be made arbitrarily small by reducing the size 
of the counting window T. An attractive feature of this representation is that it preserves 
the correspondence between the discrete time axis of the counting process {Z&} and the 
absolute "real" time axis of the underlying point process. Within the process of counts 
{Zk}, the elements and Zk +n refer to the number of counts in windows separated by 
precisely T(n — 1) sec, so that correlation in the process {Zj.} is readily associated with 
correlation in the underlying point process dN(t). 

1.2 Scaling and Power-Law Behavior in Fractal Stochastic Pro- 
cesses 

The complete characterization of a stochastic process involves a description of all possible 
joint probabilities of the various events occurring in the process. Different statistics provide 
complementary views of the process; no single statistic can in general describe a stochastic 
process completely. Fractal stochastic processes exhibit scaling in their statistics. Such 
scaling leads naturally to power-law behavior, as demonstrated in the following. Consider a 
statistic / which depends continuously on the scale x over which measurements are taken. 
Suppose changing the scale by any factor a effectively scales the statistic by some other 
factor g(a), related to the factor but independent of the original scale: 

f(ax) = g(a)f(x). (1) 

The only nontrivial solution of this scaling equation for real functions and arguments, that 
is independent of a and x is 

f(x) = bg(x) with g(x) = x c (2) 

for some constants b and c |2|, [3]]. Thus statistics with power-law forms are closely related 
with the concept of a fractal. The particular case of fixed a admits a more general solution 

& 

g(x; a) = x c cos[27r ln(x) /ln(a)] . (3) 
4 



1.3 Fractal Stochastic Point Processes (FSPPs) 



Consider, for example, a commonly encountered first-order statistic for a stochastic point 
process, the interevent interval histogram (IIH). This estimates the interevent-interval prob- 
ability density function (IIPDF) p(t) by computing the relative frequency of occurrence of 
interevent intervals as a function of interval size. This measure highlights the behavior of 
the times between adjacent events, but reveals none of the information contained in the 
relationships among these times, such as correlation between adjacent time intervals. For 
a fractal stochastic point process (FSPP) the IIPDF follows the form of Eq. (||D, so that 
p(t) ~ t c over a certain range of t, where c < — 1. 

A number of statistics may be used to describe an FSPP, and each statistic which scales 
will in general have a different scaling exponent c. Each of these exponents can be simply 
related to a more general parameter a, the fractal exponent, where the exact relation between 
these two exponents will depend upon the statistic in question. For example, the exponent 
c of the IIPDF defined above is related to the fractal exponent a by c = — (1 + a). 

Sample functions of the fractal renewal point process (FRP; see Sec. |4.1|) are true fractals; 
the expected value of their generalized dimensions (see Sec. [3.2|) assumes a nonintegral value 
between the topological dimension (zero) and the Euclidean dimension (unity) |J. 

1.4 Fractal-Rate Stochastic Point Processes (FRSPPs) 

However, the sequence of unitary events observed in many biological and physical systems 
do not exhibit power-law-distributed IIHs but nevertheless exhibit scaling in other statis- 
tics. These processes therefore have integral generalized dimensions (see Sec. p.2|) and are 
consequently not true fractals. They are nevertheless endowed with rate functions that 
are indeed either fractals or their increments: fractal Gaussian noise, fractal Brownian mo- 
tion, or closely related processes. Therefore, such point processes are more properly termed 
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fractal- rate stochastic point processes (FRSPPs). 



1.5 Relation of Fractal Exponent to Hurst Exponent 

The fractal exponent a defined above is related to the more commonly encountered Hurst 
exponent H ||. The relationship is ambiguous, however, since some authors f5j, Pi PI |10 



use the formula a = 2H + 1 for all values of a, while others |IT| use a = 2H — 1 for a < 1 
to restrict H to the range (0,1). In this paper, we avoid this confusion by considering a 
directly instead of H. 



2 APPLICATIONS OF FRACTAL AND FRACTAL- 
RATE STOCHASTIC POINT PROCESSES 

Many phenomena are readily represented by FSPPs and FRSPPs. We provide several ex- 



amples drawn from the physical [|T|, [13|, [14| and biological |15], [16], ]L7j sciences. 



2.1 Trapping Times in Amorphous Semiconductors 

A multiple trapping model has been used to show how traps that are exponentially dis- 
tributed over a large range of energies lead to a power-law decay of current in an amorphous 
semiconductor [TR [Bl EH, ETJ, EES. If a pulse of light strikes such a semiconductor, the many 



carriers excited out of their traps will be available to carry current until they are recaptured, 
which happens relatively quickly. At some point each carrier will be released from its trap 
by thermal excitation and become mobile for a time, and then be recaptured by another 
trap. For exponentially distributed energy states with identical capture cross sections, the 
electrons tend to be trapped in shallow states at first, but the probability of being caught 
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in a deep trap increases as time progresses. This leads to a current that decreases as a 
power-law function of time. 



The multiple trapping model may be usefully recast in terms of a standard fractal re- 



newal process (see Sec. 11) [23, 24]. Consider an amorphous semiconductor with localized 
states (traps) with energies which are exponentially distributed with parameter E Q between 
a minimum energy El of the order of kT, where k is Boltzmann's constant and T is the 
temperature in degrees Kelvin; and a maximum energy Eh determined by the bandgap of 
the material. For a particular trap with energy S, the corresponding mean waiting time to 
release is 

r = r exp(£/KT), (4) 

where tq is the average vibrational period of the atoms in the semiconductor. If we define 
characteristic time cutoff's A = tq exp(EL/ kT) and B = tq exp(En/ kT), and the power-law 
exponent c = k,T/Eq, then the mean waiting time r has a density which decays as a power 
law with this exponent between those two cutoffs. Each trap holds carriers for times that are 
exponentially distributed given the conditional parameter r, and averaging this exponential 
density over all possible values of r yields the unconditional trapping time density, which is 
itself approximately power law: 

p(t) « cT(c+ l)A c r {c+1 \ (5) 

for A <C t <C B. Thus each carrier will be trapped for a period that is essentially power- law 
distributed. 

Upon escaping from a trap, the carrier can conduct current for a short time until it 
is again captured by another trap. Thus each carrier executes a series of current-carrying 
jumps well described by a standard FRP. Assuming that each carrier acts independently of 
the others, the action of the carriers as a whole can be modeled as the superposition of a 
collection of such component processes, which converges to the fractal-Gaussian-noise driven 



Poisson process (FGNDP; see Sees. [4.2. 1| and f4.4|) in the limit of a large number of carriers. 
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Again, both experimental |25| and theoretical JZS, [26| results point to a power-law or fractal 
decay in the power spectral density, while several other statistics also show scaling behavior. 



2.2 Noise and Traffic in Communications Systems 



Burst noise occurs in many communications systems and is characterized by relatively brief 
noise events which cluster together, separated by relatively longer periods of quiet. Berger 
and Mandelbrot p?| , p8|j long ago showed that burst errors in communication systems are 
well modeled by a version of a fractal renewal process (see Sec. |4.1|) , and in particular that 
the interevent times were essentially independent of each other for time scales determined 
by the resolution and the duration of the observation. 



Moreover, the rate of traffic flow itself displays fractal fluctuations on a variety of high- 



speed packet-switching networks conducting different types of traffic |29|, |3(| |3T], |32], |33| • This 
has been demonstrated for time scales greater than about one sec in both the power spectral 



density and the Fano factor (see Sees. |3]J and [3.5|) . Over these time scales, the fractal-shot 



noise driven Poisson process (FSNDP; see Sees. |4.2.5| and |4.4j) has been successfully employed 



to model the traffic G3 



2.3 Current Flow in Biological Ion Channels 

Ion channels reside in biological cell membranes, permitting ions to diffuse in or out |[35| . 
These channels are usually specific to a particular ion, or group of related ions, and block 
the passage of other kinds of ions. Further, most channels have gates, and thus the channels 
may be either open or closed. In many instances, intermediate conduction states are not 



observed. Some ion channels may be modeled by a two-state Markov process [^] , with one 
state representing the open channel, and the other representing the closed channel. This 
model generates exponentially distributed dwell times in both states, which are, in fact, 
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sometimes observed. However, many ion channels exhibit independent power-law distributed 
closed times between open times of negligible duration [[37|, [3^, [39], 0, and are well described 
by a fractal renewal point process (see Sec. |Q| ) H4l[ [42 . 



2.4 Vesicular Exocytosis at the Synaptic Cleft 

Communication in the nervous system is usually mediated by the exocytosis of multiple 
vesicular packets (quanta) of neurotransmitter molecules at the synapse between two cells, 
triggered by an action potential (information-bearing signal) at the presynaptic cell [43 |. 



Even in the absence of an action potential, however, many neurons spontaneously release 



individual packets of neurotransmitter On arrival at the postsynaptic membrane, each 
neurotransmitter packet induces a miniature endplate current (MEPC); superpositions of 
MEPC-like events comprise the postsynaptic endplate currents elicited by nerve impulses 
arriving at the presynaptic cell |4l|. Analysis of the sequence of MEPC and MEPC-like 
events in a variety of preparations reveals the presence of memory over a range of times and 



frequencies [46]. The fractal-lognormal-noise driven Poisson process model (FLNDP; see 



Sees. [4.2.2| and [4.4| ) provides an excellent model for this process, for a variety of statistical 



measures 46, 47 



2.5 Action Potentials in Isolated Neuronal Preparations 

Many biological neurons carry information in the form of sequences of action potentials, 
which are localized regions of depolarization traveling down the length of an axon. Action 
potentials or neural spikes are brief and largely identical events and are well represented 
by a point process. The firing patterns of an isolated neuron and an isolated axons have 
been shown to have fractal characteristics. Musha and colleagues generated a synthetic 
spike train consisting of electrical impulses separated by independent intervals drawn from 
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the same Gaussian distribution ][48|. This nonfractal stimulus as applied to one end of 
an excised squid giant axon; the membrane voltage at the far end was considered to be 
the response spike train. The resulting power spectral density followed a l//-type decay. 
Spontaneous firing from an excised giant snail neuron yielded similar results jRJ, although 
for this preparation the data was selected to reduce fluctuations, and analyzed by interevent 
number rather than in real time. 



2.6 Auditory-Nerve-Fiber Action Potentials 



More recently, fractal behavior has been observed in the sequence of action potentials 
recorded from several in vivo vertebrate preparations. Real-time recordings were made under 
both spontaneous and driven firing conditions. Over short time scales, nonfractal stochastic 
point processes prove adequate for representing such nerve spikes, but over long time scales 



(typically greater than one sec) the fractal nature becomes evident | 50fl . Estimators of the 
rate of the process converge more slowly than for nonfractal processes, displaying fluctu- 



ations which decrease as a power-law function of the time used to estimate the rate [51|. 
With the inclusion of the refractory effects of nerve fibers, FRSPP models have been shown 
to properly characterize the statistical properties of certain potentials in peripheral primary 
auditory fibers in several species, over all time scales and for a broad variety of statistical 
measures [[40], [51], 0, |53], |54|, |56], [57], |59], |6(], [H], |62| ; only four parameters are required. 
This process possibly arises from superpositions of fractal ion-channel transitions [41, 42] or 



by fractal vesicular exocytosis in inner-ear sensory cells, as described briefly in Sees. ^3] and 
respectively. 
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2.7 Optic-Nerve-Fiber Action Potentials 



Many neurons in the mammalian visual system transmit information by means of action 
potentials as well, and FRSPPs also provide suitable models for describing the behavior of 
these neurons U63| . The gamma renewal process, which is nonfractal, has proven to be a 



useful model for some of these processes over short time scales |]64 |. However, nerve spike 
trains recorded from both cat retinal-ganglion-cell and lateral-geniculate-nucleus neurons, 
like those recorded from primary auditory neurons, exhibit fractal behavior over time scales 



greater than about 1 sec, thus necessitating the use of an FRSPP description [p3[ . Long- 
duration temporal correlation has also been observed in spike trains from cat striate-cortex 
neurons and from an insect visual interneuron IBBI . 



2.8 Central-Nervous-System Action Potentials 

Action-potential data collected from rabbit somatosensory cortical neurons and cat respi- 
ratory control neurons show experimental IIH plots which decay as power laws |67]]. In 



addition, for the data lengths considered, serial correlation coefficients among interevent in- 
tervals from these data sets do not differ from zero in a statistically significant way (However, 
differences may well emerge for longer data sets, as has been seen, for example, in auditory- 
nerve-fiber recordings ||50fl .) These two interevent-interval characteristics — independence 
and a power-law distribution — define the fractal renewal point process (see Sec. [4.1| ), which 
is an FSPP 

On the other hand, action potentials in the cat mesencephalic reticular formation (MRF) 
exhibit activity that can be modeled by an FRSPP. During REM phases of sleep, these neu- 
rons exhibit l//-type spectra []68], [7(J, as do certain hippocampal and thalamic neurons 



TTfl . The cluster Poisson point process (see Sec. fL5|) has been used successfully to model 



MRF neural activity |70|| . Furthermore, there is some evidence that these neurons have IIHs 
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whose tails decay in power- law fashion [72 1, so that a type of FSPP may prove suitable for 
describing these spike trains. 



2.9 Human Heartbeat Times 

The sequence of human heartbeats exhibits considerable variability over time, both in the 



short-term and in the long-term patterns of the times between beats |[7q, |74}|. A point process 
of the heartbeats is formed by focusing on the times between maximum contractions (R-R 
intervals). A particular FRSPP with an integrate-and-fire substrate has been constructed 



and shown to successfully describe these events [[75], |76fl . As with auditory and visual-system 
spike trains, over short time scales nonfractal point processes provide suitable models for 
the pattern of times between contractions; for times longer than roughly 10 sec, only fractal 
models suffice. Interestingly, parameters of the FRSPP used to model the data indicate 



promise for the diagnosis of various disease states of the heart |76 . 



3 ANALYSIS OF FRACTAL-RATE STOCHASTIC 
POINT PROCESSES 

As the examples described above illustrate, many phenomena are amenable to modeling 
by FRSPPs. The value of the fractal exponent a can often provide important information 
regarding the nature of an underlying process, and can also serve as a useful classification 
tool in some cases, as with human heartbeat data as mentioned in Section |2.9| . Accordingly, 
it is desirable to estimate a reliably, although this task is often confounded by a variety of 
issues H [F7|, [78 . 



Here we briefly review some of the techniques used for estimating a. By way of illus- 
tration, we apply these methods to a train of action potentials recorded from a lateral- 
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geniculate-nucleus (LGN) relay neuron in the cat visual system. There are 24285 events in 
this particular spike train, with an average interevent interval of 0.133 sec, comprising a total 
duration of 3225 sec f63fl . 



3.1 Interval Probability Density Function 

The interevent-interval histogram (IIH), an estimate of the interval probability function, 



and the generalized dimension (GD) WQ, RO, |8T|, are useful measures for fractal stochastic 



point processes, and in particular for the fractal renewal process (see Sec. |4.1|) . The IIH is 
simply the relative frequency of interevent-interval occurrence in the data set, ignoring all 
correlations among the intervals. 

Since realizations of FSPPs form true fractals, all measures considered in this paper, in- 
cluding the IIH, exhibit scaling for these processes. Thus the IIH proves useful in elucidating 
the fractal nature of FSPPs. For FRSPPs, on the other hand, the fractal structure is not 
manifested in first-order interevent-time statistics, so that the IIH does not scale over a sig- 
nificant range of interevent times. Rather, the fractal character of these processes lie in the 
dependencies among the interevent intervals, which are captured by the measures considered 
in Sees. |3~3|-|3~6- 



3.2 Generalized Dimension 

The generalized dimension is of interest for fractal stochastic point processes. If a data 
segment of length L is divided into intervals of length T, with Zk representing the number 
of points in the kth interval (see Fig. 1(c)), then the generalized dimension D q of a point 
process is defined as 



D q = lim ° VZ T , (6) 

q q-lr^o log(T) ' V ; 
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where the sum extends over all non-empty intervals. Particular cases are the capacity or box- 
counting dimension Dq, the information dimension lim^i D q , and the correlation dimension 
D2- The sum is a form of time averaging; for a stochastic point process, it is convenient 
to replace the sum by the product of L/T and the expected value of Z q . For a fractal-rate 
point process, however, analytical values of the D q will not in general equal analytical values 
of the dimensions obtained from the PG or AF, and in fact the D q will often assume integer 
values (see Sec. |L4] ). 



3.3 Coincidence Rate 

The CR measures the correlations between pairs of events with a specified time delay between 
them, regardless of intervening events, and is related to the autocorrelation function used 
with continuous processes. The CR is defined as 

x Pr {£[0, A] and £[r, r + A]} 
G{r) = hm A2 !± (7) 

where £[s, t] denotes the occurrence of at least one event of the point process in the interval 
[s, t) and t is the delay time. For an ideal fractal or fractal-rate stochastic point process 
with a fractal exponent a in the range < a < 1, the coincidence rate assumes the form 

G(t) = \6(t) + \ 2 [i + (\t I /to)*" 1 ], (8) 

where A is the mean rate of the process, S(t) denotes the Dirac delta function, and tq is a 
constant representing the fractal onset time. 

A stationary, regular point process with a CR following this form for all delay times r 
exhibits infinite mean. Further, statistics of fractal data sets collected from experiments 
exhibit scaling only over a finite range of times and frequencies, as determined by the reso- 
lution of the experimental apparatus and the duration of the experiment. Nevertheless, in 
much of the following we employ Eq. (|8|) without cutoffs since we find that the cutoffs do not 
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significantly affect the mathematical results in many cases. We employ similar ideal forms 
for other second-order measures defined later in this paper for the same reasons. 

The coincidence rate can be directly estimated from its definition. However, in prac- 
tice the CR is a noisy measure, since its definition essentially involves a double derivative. 
Furthermore, for FSPPs and FRSPPs typical of physical and biological systems, the CR 
exceeds its asymptotic value A 2 by only a small fraction at any large but practical value of r, 
so that determining the fractal exponent with this small excess presents serious difficulties. 
Therefore we do not specifically apply this measure to the LGN data, although the formal 
definition of coincidence rate plays a useful role in developing other, more reliable measures. 



3.4 Power Spectral Density 

The power spectral density (PSD) is a familiar and well-established measure for continuous- 
time processes. For point processes the PSD and the CR introduced above form a Fourier 
transform pair, much like the PSD and the autocorrelation function do for continuous-time 
processes. The PSD provides a measure of how the power in a process is concentrated in 
various frequency bands. For an ideal fractal point process, the PSD assumes the form 

S(u) = \ 2 5{u/2tt) + A [l + {u/u o y a ] , (9) 

for relevant time and frequency ranges, where A is the mean rate of events and uq is a cutoff 
radian frequency. 

The PSD of a point process can be estimated using the periodogram (PG) Sz{oj) of the 
sequence of counts, rather than from the point process itself ||82|| . This method introduces 



a bias at higher frequencies, since the fine time resolution information is lost as a result of 
the minimum count-window size. Nevertheless, since estimating the fractal exponent princi- 
pally involves lower frequencies where this bias is negligible, and employing the sequence of 
counts permits the use of vastly more efficient fast Fourier transform methods, we prefer this 
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technique. Alternate definitions of the PSD for point processes (and thus for the PG used 
to estimate them) exist; for example, a different PSD may be obtained from the real- valued 
discrete-time sequence of the interevent intervals. However, features in this PSD cannot be 
interpreted in terms of temporal frequency [[T6|j . 

Figure |2| displays the PG for the visual system LGN data calculated using the count-based 
approach. Throughout the text of this paper we employ radian frequency uj (radians per 
unit time) to simplify the analysis, while figures are plotted in common frequency f = uj/2ir 
(cycles per unit time) in accordance with common usage. For low frequencies, the PG decays 
as l/uj a , as expected for a fractal point process. Fitting a straight line (shown as dotted) to 
the doubly logarithmic plot of the PG, over the range from 0.002 Hz to 1 Hz, provides an 
estimate a ~ 0.67. 



3.5 Fano Factor 

Another measure of correlation over different time scales is provided by the Fano factor (FF), 
which is the variance of the number of events in a specified counting time T divided by the 
mean number of events in that counting time. This measure is sometimes called the index 
of dispersion of counts 0. In terms of the sequence of counts illustrated in Fig. 0(c), the 
Fano factor is simply the variance of {Z^} divided by the mean of {Z^}, i.e., 

s g[gfj__gSl ; (10) 

where E[-] represents the expectation operation. The FF may be obtained from the coinci- 
dence rate by an integral transform 

F(T) = 1 + A /V - r) \G(t) - A 2 1 dr, (11) 



AT Jo 

where A is the average rate of events of the point process, and the lower limit of the integral 
approaches zero from the positive side [pi . 
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The FF varies as a function of counting time T. The exception is the homogeneous 
Poisson point process (HPP), which is important as a benchmark in point-process theory 
just as the Gaussian is in the theory of continuous stochastic processes. For an HPP, the 
variance-to-mean ratio is always unity for any counting time T. Any deviation from unity 
in the value of F(T) therefore indicates that the point process in question is not Poisson in 
nature. An excess above unity reveals that a sequence is less ordered than an HPP, while 
values below unity signify sequences which are more ordered. For an ideal FSPP or FRSPP 
with < a < 1, the FF assumes the form 



f(t) = i + (T/nr, (12) 

where To is a fractal onset time that differs from tq in Eq. @ [see Eq. (pi)]. Therefore a 
straight-line fit to an estimate of F(T) vs. T on a doubly logarithmic plot can also be used 
to estimate the fractal exponent. However, the estimated slope of the FF saturates at unity, 



so that this measure finds its principal applicability for processes with a < 1 B 60, p3l R3 



Figure |3| shows the estimated FF curve for the same data set as shown in Fig. |2|. For 
counting times T greater than approximately 0.3 sec, this curve behaves essentially as ~ T a . 
The estimated value is a = 0.66 (dotted line), closely agreeing with the value obtained from 
the PG. Fair agreement between these two measures exists in other sensory spike-train and 
heartbeat sequence data sets |55L |63], |T6[ and in simulated FRSPPs M . 



3.6 Allan Factor 

The Allan variance, as opposed to the ordinary variance, is defined in terms of the variability 
of successive counts [^4[ [55| . In analogy with the Fano factor (FF), we define the Allan factor 
(AF) for the sequence of counts shown in Fig. [l] as 

^ - m %t ?y (13) 
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In terms of quantities defined earlier, we have 

A(T) = 2F(T) - F(2T) = l + ^= f V - r) [G(r) - G(2r)} dr. (14) 

As for the FF, the value of the Allan factor for the HPP is unity. For an ideal FRSPP with 
< a < 3, the AF also assumes a power-law- varying form 

A(T) = 1 + (T/T.r, (15) 

where T\ is another fractal onset time, different from tq and To [see Eq. (0)]; therefore a 
straight line fit of an estimate of A(T) vs. T on a doubly logarithmic plot yields yet another 
estimate of the fractal exponent. 

Figure f| shows the estimated Allan factor curve for the same data set as shown in Figs. |2] 
and H From visual inspection, this measure appears to be considerably "rougher" than the 
FF, which is typical of the data sets we have analyzed. Nevertheless it is clear that for 
counting times T greater than approximately 0.3 sec, its behavior can also be approximated 
as ~ T a . To estimate the value of a, a straight line fit to the doubly logarithmic plot of the 
estimate of A(T) vs. counting time T was provided. The value of a = 0.65 obtained agrees 
well with the values calculated using the PG and FF. Excellent agreement between the AF 



and the PG exists in other sensory spike-train and heartbeat sequence data sets [60. 53. 7q, B6 



and simulations |63 . Furthermore, the AF appears to yield superior estimates of a than the 



FF in all cases pH, |6B], |76|, [56|. In particular, instead of saturating at a power-law exponent 



of unity, estimates of a obtained from the AF can range up to a value of three |60|| . Thus 
we do not employ the FF in the remaining sections of this paper. 

The Allan variance, E [(Zk+i — Zk) 2 } may be recast as the variance of the integral of the 
point process under study multiplied by the following function: 

' -1 for -T < t < 0, 

^Haar(t) = | +1 for < t < T, (16) 

otherwise. 

V 

Equation (0) defines a scaled wavelet function, specifically the Haar wavelet. This can 
be generalized to any admissible wavelet ip(t); when suitably normalized, the result is a 



wavelet Allan factor (WAF) |86], This generalization enables an increased range of 

fractal exponents to be estimated, at the cost of reducing the range of times over which 
the WAF varies as T a . In particular, for a particular wavelet with regularity (number of 
vanishing moments) R, fractal exponents a < 2R + 1 can be reliably estimated |86| , |37|j . For 
the Haar basis, R — 1, while all other wavelet bases have R > 1. Thus the WAF employing 
bases other than the Haar should prove useful in connection with FRSPPs for which a > 3; 



for processes with a < 3, however, the AF appears to be the best choice [86, 87 



3.7 Relationships Among the Measures 

For an ideal FSPP or FRSPP with < a < 1, any one of Eqs. (|), (|), (03), and (0) implies 
the other three, with P, [88[ 

\r^ a T« = a(a + l)/2 

c^T a = cos(7ra/2)r(a + 2) (17) 
T a T - a = 2-2°. 

For larger values of a, the FF and CR do not scale, although the PSD and AF do; thus 
Eqs. (|9]), and (|T5| ) imply each other for these FRSPPs. Therefore, over the range of a for 
which the AF exhibits scaling [ S8| 



cos(a7r/2)r(2 + a)/{2 - 2 a ) for < a < 1 
u*T? = { vr/ln(4) for a = 1 (18) 

-cos(a7r/2)r(2 + a)/(2 Q -2) for 1 < a< 3. 

In principle, any of the these statistics could be solved for a when it is known that a < 1, 
although best results obtain from the PG and the AF. 

Since the GD of a fractal cannot exceed the Euclidean dimension, which assumes a value 
of unity for one- dimensional point processes, FSPPs cannot exist for a > 1; only FRSPPs 
are possible at these values of the fractal exponent. 
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4 SYNTHESIS OF FRACTAL AND FRACTAL- RATE 
STOCHASTIC POINT PROCESSES 



In previous work we defined an FSPP and several FRSPPs and derived a number of their 



statistics [EL E3L 891. We include brief summaries of these here, as well as of the clustered 



Poisson point-process model of Griineis and colleagues K|. We also introduce two new 
fractal-rate processes (fractal lognormal noise and fractal exponential noise) and define two 
new methods for generating a point processes from a rate function: integrate-and-fire and 
controlled variability integrate-and-fire. These prove useful in isolating the effects of fractal 
components from those of nonfractal noise introduced through the mechanics of generating 
the point process. For all of the processes considered, the scaling relation in Eq. (0) holds 
for a number of statistical measures. 



4.1 Fractal Renewal Point Process (FRP) 

The one-dimensional homogeneous Poisson point process (HPP) is perhaps the simplest 



nonfractal stochastic point process [|9T]. The HPP is characterized by a single constant 
quantity, its rate, which is the number of events expected to occur in a unit interval. A 
fundamental property of the HPP is that it is memoryless; given this rate, knowledge of 
the entire history and future of a given realization of a HPP yields no information about 
the behavior of the process at the present. The HPP belongs to the class of renewal point 
processes; times between events are independent and identically distributed. The HPP is the 
particular renewal point process for which this distribution follows a decaying exponential 
form. 

We now turn to a renewal process that is fractal: the standard fractal renewal process 



(FRP) [p_3|, EJ], Ej], E_8], f|T|. In the standard FRP, the times between adjacent events are 



independent random variables drawn from the same fractal probability distribution. In 
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particular, the interevent-interval probability density function p(t) decays essentially as a 
power law function 

= { k /t a+1 for A <t < B 
1 otherwise, 

where a is the fractal exponent, A and B are cutoff parameters, and k is a normalization 
constant. The FRP exhibits fractal behavior over time scales lying between A and B. This 
process is fully fractal for < a < 1: the power spectral density, coincidence rate, Fano 
factor, Allan factor, and even the interevent-time survivor function all exhibit scaling as in 
Eq. (0) with the same power-law exponent a or one simply related to it. 

Further, for this process the capacity or box-counting dimension D assumes the value a 
23fl ; since the FRP is ergodic, the generalized fractal dimension D q becomes independent of 



the index q [81], so D q = D = a for all q, and all fractal dimensions coincide. 



A different (nonrenewal) point process results from the superposition of a number K 
of independent FRPs; however, for this combined process as K becomes large, and indeed 
for any FRSPP, the interevent-time probability density function no longer scales, and the 
generalized dimensions D q no longer equal a, although the PG and AF (and the FF and 
CR, for < a < 1) retain their scaling behavior. As the number of such processes increases, 
and for certain ranges of parameters, this superposition ultimately converges to the fractal 
Gaussian-noise driven Poisson process [^, |1J (see also Sees. [4.2.1| and f4.4| ). 



The standard FRP described above is a point process, consisting of a set of points or 
marks on the time axis as shown in Fig. |5|(a); however, it may be recast as a real- valued 
process which alternates between two values, for example zero and unity. This alternating 
FRP would then start at a value of zero (for example), and then switch to a value of unity at 
a time corresponding to the first event in the standard FRP. At the second such event in the 
standard FRP, the alternating FRP would switch back to zero, and would proceed to switch 
back and forth at every successive event of the standard FRP. Thus the alternating FRP 
is a Bernoulli process, with times between transitions given by the same interevent-interval 
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probability density as in the standard FRP, as portrayed in Fig. ||(b). 



4.2 Integrate-and-Fire Fractal-Rate Stochastic Point Processes 

Many point processes derive from a continuous-time function \(t) which serves as the stochas- 
tically varying rate of the point process. These are known as fractal-rate stochastic processes 
(FRSPPs). We confine our discussion to rate processes (and therefore FRSPPs) which are 
ergodic. Perhaps the simplest means for generating a point process from a rate is the 
integrate-and-fire (IF) method ||76|1 . In this model, the rate function is integrated until it 
reaches a fixed threshold 8, whereupon a point event is generated and the integrator is reset 
to zero. Thus the occurrence time of the (k + l)st event is implicitly obtained from the first 
occurrence of 

[ tk+1 X(t) dt = 6. (20) 
Jt k 

With such a direct conversion from the rate process to the point process, any measure 
applied to these point processes will return results closely related to the fractal structure of 
the underlying rate process A(t). In particular, over small frequencies, the theoretical PSDs 
of the rate process and of the resulting point process coincide. 

We now turn to methods for generating several different kinds of fractal rate functions. 
More complex methods for point-process generation, as delineated in Sees. |4.jj| and |4.4j , can 
also be applied to these same rate functions. 



4.2.1 Fractal Gaussian Noise (FGN) 



Gaussian processes are ubiquitous in nature and are mathematically tractable. The values 
of these processes at any number of fixed times form a jointly Gaussian random vector; 
this property, the mean, and the spectrum completely define the process. For use as a rate 
process, we choose a stationary process with a mean equal to the expected rate of events, 
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and a spectrum of the form of Eq. (jfty with cutoffs. FGN properly applies only for a < 1; 
the range 1 < a < 3 is generated by fractal Brownian motion ||. However, in the interest 
of simplicity, we employ the term fractal Gaussian noise for all values of a. A number of 
methods exist for generating FGN || [92|, [53], |94|]: typically the result is a sampled version 
of FGN with equal spacing between samples. Interpolation between these samples is usually 
achieved by selecting the most recently defined value. 

With a rate process described by FGN serving as the input to an IF process, the fractal- 
Gaussian-noise driven integrate-and-fire process (FGNIF) results. Generally, we require that 
the mean of the rate process be much larger than its standard deviation, so that the times 
when the rate is negative (during which no events may be generated) remain small. The 
FGNIF has successfully been used to model the human heartbeat |75], |76| (see Sec. |2li| ). 

All of the rate processes considered in the following subsections converge to FGN under 
appropriate circumstances, and thus the point processes they yield through the IF construct 
converge to the FGNIF. 



4.2.2 Fractal Lognormal Noise (FLN) 



A related process results from passing FGN through a memoryless exponential transform. 
Since the exponential of a Gaussian is lognormal, we call this process fractal lognormal noise 
(FLN). If X(t) denotes an FGN process with mean E[X], variance Var[X], and autoco- 
variance function Kx(t), then X(t) = exp[X(£)] is an FLN process with moments E[A n ] = 
exp{nE[A] + n 2 Var[X]/2} and autocovariance function K\{r) = E 2 [X]{exp[Kx(T)] — 1}- 

When this rate serves as the input to a doubly stochastic Poisson point process (see 



Sec. |4~4| ), the result provides an excellent model for vesicular exocytosis (see Sec. |2~4] ). By 
virtue of the exponential transform, the autocorrelation functions of the Gaussian process 
X(t) and the lognormal process X(t) differ; thus their spectra differ. In particular, the PSD 
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of the resulting point process does not follow the exact form of Eq. (|9]), although for small 
values of Var[X], as encountered in modeling exocytosis (for which Var[X] < 0.6) |46|, f47f , 
the experimental PG closely resembles this ideal form. In the limit of very small Var[X], the 
exponential operation approaches a linear transform, and the rate process becomes FGN. 

4.2.3 Fractal Exponential Noise (FEN) 

Other nonlinear transforms may be applied to FGN, yielding other fractal rate processes; we 
consider one which generates an exponentially distributed process. If Xi(t) and X 2 (t) denote 
two independent and identically distributed FGN processes with zero mean, variance Var[X], 
and autocovariance function Kx{t), then X(t) = Xf(t) + X|(t) is a fractal exponential 
noise (FEN) process with moments E[A n ] = 2 n n\ (Var[X]) n and autocovariance function 
K\(r) = 4K x (t). The rate X(t) turns out to be exponentially distributed. If Kx{t) scales 
as in Eq. with an exponent ax in the range 1/2 < ax < 1, then so will K\(t), but with 
an exponent a\ = 2a x — 1. 

This process may prove useful in the study of certain kinds of thermal light ||95|| . Such light 
has an electric field which comprises two independent components, with Gaussian amplitude 
distributions. The light intensity is the sum of the squares of the two fields, and therefore 
has an exponential amplitude distribution. 



4.2.4 Fractal Binomial Noise (FBN) 

A number of independent, identical alternating fractal renewal processes (see Sec. [4.ip may be 



added together, yielding a binomial process with the same fractal exponent as each individual 
alternating FRP p3| . This binomial process can serve as a rate process for an integrate- 



and-fire process; the fractal-binomial-noise driven integrate-and-fire (FBNIF) results. It is 
schematized in Fig. ^ As the number of constituent processes increases, the FBN rate 
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converges to FGN (see Sec. 



4.2.1 ) with the same fractal exponent. 



4.2.5 Fractal Shot Noise (FSN) 



Though the HPP itself is not fractal, linearly filtered versions of it, denoted shot noise, can 
exhibit fractal characteristics. In particular, if the impulse response function of the filter 



has a decaying power-law (fractal) form, the result is fractal shot noise (FSN) |96|, £)7], eE. 



If FSN serves as the rate for an IF process, the fractal-shot-noise driven integrate-and-fire 
process (FSNIF) results; Fig. [7] schematically illustrates the FSNDP as a two-stage process. 
The first stage produces a HPP with constant rate //. Its output M(t) becomes the input to 
a linear filter with a power-law decaying impulse response function 

h(t) = { k ^ a/2 for A <t<B 
I otherwise, 

where a is a fractal exponent between zero and two, A and B are cutoff parameters, and k 
is a normalization constant. This filter produces fractal shot noise I(t) at its output, which 
then becomes the time-varying rate for the last stage, an integrate-and-fire process. The 
resulting point process N(t) reflects the variations of the fractal-shot-noise driving process. 
Under suitable conditions, FSN converges to FGN as provided by the central limit theorem 
1; the result is then the FGNIF. 



4.3 Controlled- Variability Integrate-and-Fire Fractal-Rate 
Stochastic Point Processes 

FRSPPs based on an integrate-and-fire substrate have only one source of randomness, which 
is the rate process A(i). Those based on a Poisson-process substrate (see Sec. |4.4|) have a 
second source which depends explicitly on the rate process. We now generalize the family of 
FRSPPs to include those which have a second source of randomness that may be specified 
independently of the first. This new process differs from the simple IF process based on 
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X(t), by the imposition of a specific form of jitter on the interevent times. After generating 
the time of the (n + l)st event t n+ \ in accordance with Eq. (P0|) , the (n + l)st interevent 
time t n+ i — t n is multiplied by a Gaussian-distributed random variable with unit mean and 
variance a 2 . Thus t n+ i is replaced by 

t n+1 + a (t n+1 -t n )Af (0,1), (22) 

where A/"(0, 1) is a zero-mean, unity-variance Gaussian random variable. The result is the 
jittered-integrate-and-fire (JIF) family of processes. Employing the fractal-rate functions de- 
lineated in Sec. |4.2| yields the fractal-Gaussian-noise driven jittered-integrate-and-fire process 
(FGNJIF) (Sec. [4.2. 1|) , the fractal-lognormal-noise driven jittered-integrate-and-fire process 
(FLNJIF) (Sec. |4.2!2| ), the fractal-exponential-noise driven jittered-integrate-and-fire process 



(FEN JIF) (Sec. [4.2.3Q , the fractal-binomial-noise driven jittered-integrate-and-fire process 



(FBNJIF) (Sec. 4.2.4j ) , and the fractal-shot-noise driven jittered-integrate-and-fire process 



(FSNJIF) (Sec. 



In these point processes, the standard deviation o of the Gaussian jitter is a free parameter 
that controls the strength of the second source of randomness. Two limiting situations exist. 
In the limit o —>■ the second source of randomness is absent and the JIF processes reduce 



to the simple IF process considered in Sec. 4.2. The opposite limit, a — ► oo leads to a 



homogeneous Poisson process; none of the fractal behavior in the rate process appears in 
the point process, as if the rate were constant. Between these two limits, as a increases 
from zero, the fractal onset times of the resulting point processes increase, and the fractal 
characteristics of the point processes are progressively lost, first at higher frequencies (shorter 
times) and subsequently at lower frequencies (longer times). 

Finally, we note that another generalization of the integrate-and-fire formalism is possible 
||76|| . The threshold for firing 9 need not remain constant, but can itself be a random variable 
or stochastic process. Randomness imparted to the threshold can then be used to represent 
variability in other elements of a system, such as the amount of neurotransmitter released per 
exocytic event or the duration of ion channel openings. The homogeneous Poisson process 
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(HPP) with rate Ao would then form a special case within this construct, generated when the 
rate is fixed and constant at A(i) = Ao, and 9 is an independent, exponentially distributed 
unit mean random variable associated with each interevent interval. If the rate process X(t) 
is also random, a doubly stochastic Poisson point process results instead of the HPP. 



4.4 Fractal-Rate Doubly Stochastic Poisson Point Process 
(FRDSPP) 

We have seen that an integrate- and- fire process (with or without jitter) converts a fractal 
rate process into an FRSPP. A Poisson process substrate may be employed instead of the 
IF substrate, in which case the outcome is the family of doubly stochastic Poisson point 



processes (DSPPs) [9J|. DSPPs display two forms of randomness: one associated with the 
stochastically varying rate, and the other associated with the underlying Poisson nature of 
the process even if the rate were fixed and constant. With fractal rate processes, a fractal- 
rate DSPP (FRDSPP) results 0. Again, simple relations exist between measures of the rate 
X(t) and the point process dN(t). In particular, for the theoretical PSDs we have 

S N (u) = E[\] + S x (u), (23) 

where Sn(u) is the PSD of the point process and S\(u) is the PSD of the rate. 

As with the fractal IF processes considered above, the amount of randomness introduced 
by the second source in fractal DSPPs (the Poisson transform) may also be adjusted. How- 
ever, in this case the amount introduced depends explicitly on the rate A(i), rather than 
being independent of it as with the JIF family of FRSPPs. For example, for a particular 
rate function A(t) suppose that its integral between zero and t assumes the value A: 

' \{u) du = A. (24) 



o 



Then given A, the number of events N(t) generated in the interval (0, t] has a Poisson 
distribution, and in particular the mean and variance of N(t) both assume a value of A. 
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Therefore, if X(u) is multiplied by a constant value, then for a fixed integration time t, the 
mean and variance of N(t) will change by the same factor. As this constant value increases, 
the coefficient of variation (the standard deviation divided by the mean) of N(t) decreases, 
and the contribution of the Poisson transform to the total randomness of the point process 
decreases, ultimately becoming negligible as the multiplication constant increases without 
limit. However, the DSPP formalism is less flexible than the JIF formalism for introducing 
randomness, since in the former case it explicitly depends on the mean number of events, 
while in the latter case the variance may be independently adjusted with the parameter a. 

Again, a variety of fractal point processes exist in the FRDSPP class, depending on 
the rate process chosen. Examples include the fractal-Gaussian-noise driven Poisson point 
process (FGNDP) (Sec. [4.2.1| , ||89||), the fractal-lognormal- noise driven Poisson point pro- 
cess (FLNDP) (Sec. [4.2.2| , [[46| , f47j), the fractal-exponential- noise driven Poisson point pro- 



cess (FENDP) (Sec. |4.2.3[ ) , the fractal-binomial-noise driven Poisson point process (FB- 
NDP) (Sec. P~2~4j HI), and the fractal-shot -noise driven Poisson point process (FSNDP) 
(Sec. |4.2.5| , ||89|| ). For the FGNDP a slight modification is necessary since a negative rate 
has no meaning for a Poisson point process. Therefore, rather than using a Gaussian process 
directly we instead set the rate to zero for negative values of this process. We require that 
the mean rate be much larger than its standard deviation, so that the effect of this nonlinear 
limiting can essentially be ignored. Again, all forms of FRSPPs with positive finite mean 
rates converge to the FGNDP under appropriate conditions 0. 



4.5 Cluster Poisson Point Process 

Other important formulations for FRSPPs exist. Griineis and colleagues defined a clustered 
Poisson point process in which each member of a sequence of primary events from an HPP 
gives rise to a train of secondary events (as in the FSNDP), but where the events in the 
train form a segment of a renewal process (usually also an HPP, but often with a different 
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rate), with a fractal (power-law distributed) duration ||100| , |90|j . The resulting process indeed 



exhibits power-law scaling in the same statistics as other FRSPPs [|101|1 . The cluster Poisson 
point process is therefore a Bartlett-Lewis-type process, whereas the FSNDP is a process of 
the Neyman-Scott type This process has been used successfully to model mesencephalic 
reticular formation neural spike trains []70[ (see Sec. |2.8| ). 



5 ESTIMATION OF THE ALLAN FACTOR AND 
POWER SPECTRAL DENSITY: THEORY 

Given a segment of an FRSPP, we seek an estimate a of the true fractal exponent a of the 
entire process from which the segment was derived. Several effects contribute to estimation 
error for finite-length segments, regardless of the method used. The fractal exponent provides 
a measure of the relative strengths of fluctuations over various time scales; for an FRSPP 
with a relatively large fractal exponent, for example, relatively more energy is concentrated in 
longer time scale fluctuations than for an FRSPP with a smaller fractal exponent. Variance 
stems from the inherent randomness of the strengths of these fluctuations. A collection of 
finite realizations of an FRSPP with the same parameters will exhibit fractal fluctuations of 
varying strengths, leading to a distribution of estimated fractal exponents. 

Bias appears to arise from cutoffs in the time (frequency) domain which give rise to 
oscillations in the frequency (time) domain, confounding the pure power-law behavior of 
the fractal PG (AF). In addition, the physical limitations of the measurement process itself 
impose practical limits on the range of time scales available. Although algorithms exist 
for accurately compensating for these cutoffs, they presuppose a detailed knowledge of the 
process a priori, which is not in the spirit of estimating an unknown signal. Consequently, 
in this paper we do not attempt to compensate for these cutoffs in this manner. In previous 
work we determined the theoretical expected Fano factor values for an FRSPP with finite 
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duration, and the corresponding expected bias for the corresponding fractal exponent. Bias 
for the power spectral density was also considered, but not for a finite data length. 

The AF and the PG highlight the fractal behavior of FRSPPs particularly well, and thus 
prove to be most useful for estimating fractal exponents, as indicated in Sec. |3|. We proceed 
to investigate the bias in these measures, and the effects of this bias on fractal-exponent 
estimation. The variance of PG-based estimators was examined previously ||; that for the 
AF appears not to be readily amenable to analytical study. 



5.1 Effects of Finite Data Length on the AF Estimate 

Unlike the Fano factor, estimates of the AF do not suffer from bias due to finite data length; 
thus fractal exponent estimates based on the AF do not either. Given a particular data set, 
the estimated AF A(T) at a particular counting time T is given by 

N-2 I N-l 

A(T) — (N — l)- 1 J2 (Z k+ i - Z k f / 2N- 1 J2 Z k , (25) 

k=0 I k=0 

where N is the number of samples, {Zk} represents the sequence of counts as shown in 
Fig. |l](c), and the functional dependence of Zk upon T is suppressed for notational clarity. 

This estimate of the AF is simply the estimate of the Allan variance divided by twice 
the estimate of the mean; however, computing the expected value of this estimate is not 
straightforward. We therefore employ the true mean rather than its estimate in Eq. (f2~5|). 
This does not appreciably affect the result, since the error so introduced remains a constant 
factor for all counting times, and so cancels in power-law slope calculations where logarithms 
are used. In any case, the estimate of the Allan variance exhibits far larger variations than 
the estimate of the mean, so the fluctuations in the AF estimate are dominated by the 
former. Therefore, the expected value of the AF estimate becomes 



E 



A(T)] = E{E[(Z k -Z k+1 ) 2 ]/2E[Z, 
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E 
E 



(Z - Z 1 ) 



/2E[Z k ] 
2E[Z ] 



(AT)- X E \ZI - Z Q Z 1 



(26) 



which is independent of the number of samples, N, and therefore of the duration of the 
recording L. In the limit L — > oo, the expected estimated AF for an FRSPP with < a < 3 
follows the form of Eq. (|15|) precisely, by definition of an ergodic process. Since Eq. ( f26|) 
does not depend on L, the quantity E A(T) must assume a constant value independent 
of L; thus the AF does not exhibit bias caused by finite data duration. Finally, since the 
expected value of the AF estimate is unbiased, the expected value of the estimate of a itself 
is expected to have negligible bias. This is a simple and important result. 

Estimates of the ordinary variance and Fano factor, in contrast, do depend on the dura- 
tion. In this case we have 



Var[Z] 



N-l 

(N-i)-^(z k -nz]r 

k=0 

N-l / N-l \ 

(tf-ir'E [Zk-N-^zA 

k=0 \ 1=0 / 



(N-l) 



-i 



N-l 



N-l N-l 



N-l N-l N-l 



E zl - £ £ z k z x + iv- 2 E E E w« 

.k=0 k=0 1=0 k=0 1=0 rn=0 

N-l N-l N-l 

(N - I)' 1 E Z 2 k - N-\N - l)- 1 £ £ Z k Z, 

k=0 k=0 1=0 

N-l N-l 

N-^Z 2 k -N-\N-l)^j:j:ZkZ, 

k=0 k l^k 
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These cross terms, which do depend on the number of samples, lead to an estimated Fano 
factor with a confounding linear term 



E 



F(T) w 1 + (T/T ) a — T [TqL 



-to. t 1— a 



(28) 



for < a < 1 |2J. The last term on the right-hand-side of Eq. leads to bias in estimating 
the fractal exponent; for this and other reasons, we do not employ the FF in fractal-exponent 
estimation. 
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5.2 Effects of Finite Data Length on the PSD Estimate 



Computation of the periodogram (estimated PSD) proves tractable only when obtained 
from the series of counts {Zk}, rather than from the entire point process N(t). This, and 
more importantly the finite length of the data set, introduce a bias in the estimated fractal 
exponent as shown below. 



We begin by obtaining the discrete-time Fourier transform of the series of counts {Z k }, 
where < k < M: 

M-1 

-j2Trkn/M 



M-1 

Z n = E 

k=0 



(29) 



The periodogram then becomes 



S 7 (n) = M- 1 \Z U 



n 

M-1 



M" 1 ]T Z k Z m e^ k ~ m ^ M , 

k=0 m=0 
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with an expected value 
E 



M-1 M-1 

S z (n)\ = M- 1 J2 E e^ k - m ^ M E[Z k Z m ] 

k=0 m=0 

M-1 M-1 I" -t r T 

= M~ l V V e j27T{k - m)n/M E / / dN(s + kT)dN(t + mT) 

Js=oJt=0 



k=0 m=0 
M-1 M-1 



T r T 



= M ~ X E E e f2 ^ k - m)n/M I I G N [s-t+{k- m)T\ ds dt 
k=0 m=0 Js=0 Jt=0 

M ~ lM - 1 , rT r2T-\u\ dll dl) 

= ikT 1 ]T £ e^ k ~ m ^' M / / G N [u + (k - m)T\ — 

k=0 m=0 Ju=-TJv=\u\ 2 

M-1 M-1 r T 



= M' 1 Y, E e j2n(k ~ m)n/M f _ (T - |u|) G N [u + (k - m)T] du 

M-1 M-1 r T 

= M- 1 E e j2 ^ k ~ m)n/M / (T — \u\) 

k=0 m=0 Ju=-T 

x r S N {u)e^ k - m ^ d ^du 

Jlu=-oo 2n 

= (27VM)- 1 / S 



N{UJ 



M-1 
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jk(2wn/M+u)T) 



u=-T 



(T — \u\) e ju)U dudu 
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(2nM)~ l 

T 

nM 



_ , , sm 2 (7m + MuT/2) 4sin 2 (u/772) , 
o 7v {ujj— — — — — — n uu; 



sin i (7rn/M + c«;T/2) a/ 



S N (2x/T) 



sin 2 (Ma;) sin 2 (a;) 
x 2 sin 2 (x + nn/M) 



dx. 
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For a general fractal stochastic point process, where the PSD follows the form of Eq. 
we therefore have 
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x 2 sin 2 (a; + nn/M) 



dx. 



(32) 



Focusing on the smaller values of n useful in estimating fractal exponents permits the use 
of two approximations in Eq. fl32|). For the values of a used in this paper (0 < a < 2), the 
integrand of Eq. (^) will only be significant near x = —irn/M, yielding 

AT r°° r. , , „.i MnS(x + nn/ M) sin 2 (x) 



E 



S z n) 



TlM J-oc 

AT ' 



1 + (Lu T/2) a \x\ 
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1 + (2irn/u MT)- a 
1 + (2im/u MT)- a 



sin 2 (nn/M) 
(nn/M) 2 



(33) 



which is of the same form as Eq. (Q). Improvement of this estimation procedure appears to 
require numerical integration of Eq. ([32]), which proves nontrivial since the integrand exhibits 
oscillations with a small period n/M. Fortuitously, for the parameter values employed in 
this paper the integrand appears peaked near x = —nn/M, so that not too many (« 200) 
oscillations need be included in the calculations. Indeed, numerical results employing this 
method, and with uiq — > for which Eq. ([32]) is known to assume the simple value AT, agree 
within 0.1%. These results form an extension of earlier approaches [Q, Q which ignored the 
effects of imposing periodic boundary conditions on the Fourier transform, and of binning 
the events. Numerical integration of Eq. (^) followed by a least-squares fit on a doubly 
logarithmic plot leads to results for the expected bias of the PSD-based estimate of the 
fractal exponent, shown in Table |[ 



Other methods exist for compensating for finite data length in power spectral estimation, 



33 



such as maximizing the entropy subject to the correlation function over the available range 
of times (see, e.g., ||102|| ). 



6 ESTIMATION OF THE ALLAN FACTOR AND 
POWER SPECTRAL DENSITY: SIMULATION 



Having considered various possible theoretical sources of error in simulating FRSPPs with a 
desired fractal exponent, we now proceed to investigate their effects upon fractal exponent 
estimates based on the power spectral density and Allan factor measures. To this end 
we employ simulations of three of the FRSPPs outlined above: FGN driving a deterministic 
integrate-and-fire process (FGNIF), FGN driving an IF process followed by Gaussian jitter of 
the interevent intervals (FGNJIF), and FGN driving a Poisson process (FGNDP). We choose 
these three processes from the collection presented in Sec. |] for a number of reasons. First, 
the FGNIF, FGNJIF, and FGNDP derive from a common continuous-time rate process, 
fractal Gaussian noise (FGN), and thus differences among these point processes must lie in 
the point-process generation mechanism rather than in the fractal rate processes. Second, 
FGN admits fractal exponent of any value, while the fractal exponents cannot exceed two 
for fractal binomial noise and fractal shot noise. Third, of all the FRSPPs examined in 
Sec. |], those based on the FGN appear to suffer the least from the effects of cutoffs, so 
that expected values of the PG and the AF most closely follow the pure power-law forms of 
Eqs. (H) and ([To]) respectively. 

To generate the rate function with 1/ f a spectral properties, we computed discrete-time 
approximations to fractal Gaussian noise (FGN) by forming a conjugate-symmetric sequence 
X[k] of M samples which satisfies 



where p = E[A(t)] is the desired average value of the FGN process, and c is a constant 




for k = 

for 1 < k < M/2, 



(34) 
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determining the relative strength of the FGN. The phases of X[k] were independent and 
uniformly distributed in the interval [0, 2ir) for 1 < k < M/2, while both X[0] and X[M/2] 
were strictly real. A discrete-time sequence x[n] was obtained by taking the inverse discrete 
Fourier transform of X[k]. The resulting FGN process is periodic with period M, so to reduce 
the effects of this periodicity we employ only the first half of the sequence. Specifically, 
we used an original sequence x[n] of length M = 2 16 = 65 536 points and kept the first 
2 15 = 32 768 samples as x[n]. Further, since rate functions must not attain negative values, 
we chose the mean value of the FGN, p, sufficiently large compared to the relative strength 
of the FGN c to ensure that x[n] remained positive for all 2 15 different values of x[n] in all 
simulations. In particular for a < 1 this was adjusted by choosing the FF intercept time T 
to be 10 times the average interevent time, while for a > 1 the PSD fractal onset frequency 
ujq was chosen to be 0.0005 times the average rate p. 

The resultant sequence x[n] was taken to represent equally spaced samples of a continu- 
ous-time process; to simplify subsequent analysis, we chose the duration of each sample to be 
one sec, without loss of generality. Varying the mean event production rate p while keeping 
all other parameters fixed (or in fixed relation to p, as with the FF intercept time and PSD 
corner frequency) resulted in point processes which contained varying expected numbers of 
events. In particular, we employed values of p = 5, 10, 20, 40 and 80 events per sample of 
x[n], leading to a range of roughly 164 000 to 2 620 000 events in the resulting point process. 

The following procedures were employed to compute the estimates described in Sec. [| 
For the PG, the absolute event times were quantized into 2 16 bins, which therefore formed 
a rate estimate of the process with counting windows having a duration of 0.5 sec. Then a 
discrete Fourier transform was performed, followed by replacing the Fourier components with 
their square magnitudes. Finally, a least-squares fit was performed on the logarithm of these 
quantities vs. the logarithm of the frequencies / for various selected ranges of frequencies. 
These slopes form PG-based estimates of the fractal exponent. 

Estimates of the AF were calculated over a range of counting times T in a geometric 
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sequence with ten counting times per decade. The slopes of the AF estimate were again 
obtained from linear least-squares fits to the logarithm of the AF values vs. the logarithm of 
the size of the counting time T; these slopes became the AF-based estimates of the fractal 



exponent. As indicated in Sec. |5J], finite data lengths do not impart extraneous bias to the 
AF-based estimates. 

For each type of FRSPP and value of p employed, 100 simulations were performed for 
each of three different design values of the fractal exponent: a = 0.2, 0.8, and 1.5; AF 
and PG curves were computed for each simulation. Estimates of the fractal exponent were 
obtained by two methods for each of these two measures. In the first, all 100 AF or PG curves 
were averaged together, after which the doubly logarithmic slope of this average curve was 
computed, yielding a single estimate of the fractal exponent. The second method consisted 
of first obtaining the doubly logarithmic slopes of the individual AF or PG curves, and then 
averaging these slopes. This second method yields the mean fractal exponent as well as the 
standard deviation (and other possible statistics of the fractal exponent). 



As shown in Sec. ^7|, the FGN-driven integrate-and-fire process (FGNIF) is obtained by 
identifying the FGN sequence x[n] as samples of a rate function X(t), taken to be constant 
over each sample of x[n] (1 sec). This rate function is integrated until the result attains 
a value of unity, whereupon an output point is generated and the integrator is reset to 
zero. The generation of the continuous rate function from the piecewise- const ant FGN does 
not significantly change the observed fractal structure of the resulting point process. For 
frequencies less than the inverse of the sampling time T (1 Hz in this example), the PSDs of 
the piecewise-constant and exact versions of FGN differ by a factor of sinc 2 (u;T/2); fractal 
behavior depends on low frequencies (u <C 1/T) where the above factor is essentially unity. 
Since the FGNIF process has just one source of randomness, namely the underlying FGN, 
the estimators for the fractal exponent should therefore display only the (fractal) behavior 
of the FGN itself, along with finite-size effects. Bias or variance due to a second source of 
randomness, such as the Poisson transform in the FGNDP, will not be present. Thus for 
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our purposes the FGNIF process serves as a benchmark for the accuracy of fractal exponent 
estimators, with a bias and variance depending only on finite data size effects and not on 
the rate p. 

The results for the FGNIF process with a rate p = 40 and three different fractal exponents 
(ct=0.2, 0.8 and 1.5) are summarized for the AF in Table [l] and for the PG in Table ^|. 



Range (sec.) 


Method 


a = 0.2 


a = 0.8 


a = 1.5 


62.5 - 625 


Fit of Avg. 
Avg. of Fits 


0.199 

0.194 ±0.074 


0.799 

0.795 ± 0.072 


1.499 

1.495 ±0.067 


125 - 1250 


Fit of Avg. 
Avg. of Fits 


0.211 

0.199 ±0.119 


0.807 

0.807 ±0.114 


1.499 

1.490 ±0.105 


250 - 2500 


Fit of Avg. 
Avg. of Fits 


0.219 

0.184 ±0.165 


0.804 

0.804 ±0.153 


1.490 

1.472 ±0.138 


25 - 2500 


Fit of Avg. 
Avg. of Fits 


0.204 

0.192 ±0.056 


0.801 

0.792 ±0.057 


1.495 

1.487 ±0.059 



Table 1: AF-based fractal exponent estimates a for simulated FGNIF processes with a 
rate p = 40, for different time ranges. The governing rate processes have theoretical fractal 
exponents of a = 0.2, 0.8, and 1.5. The estimated fractal exponents were obtained from 
straight-line fits on doubly logarithmic coordinates to an average curve of 100 independent 
simulations (Fit of Avg.) and from averages of the slopes of the individual curves (Avg. 
of Fits). The deviations are the standard deviations obtained from the second averaging 
procedure. 

These Tables present estimated values of the fractal exponent a obtained from averages over 
100 independent simulations, for a variety of time and frequency ranges. The exponents were 
obtained by the two methods delineated above, and are labeled "Fit of Avg." and 'Avg. of 
Fits" , respectively; theoretical results for the PSD are labeled "Theory" . For an experimental 
AF or PG following the form of Eq. (|T^) or ([5J), respectively, the range of times or frequencies 
over which the fractal exponent is estimated will not significantly affect the returned value, 
as long as the range lies in the scaling region ujq 1 <^T<^L (or 1/L<tiu)<^ u ), respectively. 
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Range (Hz) 


Method 


a = 0.2 


a = 0.8 


a = 1.5 


0.00025 - 0.0025 


Fit of Avg. 


0.200 


0. 


335 


1.628 




Avg. of Fits 


0.201 ±0.140 


0. 


836 ±0.137 


1.623 ±0.118 




Theory 


0.200 ±0.116 


0. 


307 ±0.116 


1.545 ±0.116 


0.0005 -0.005 


Fit of Avg. 


0.196 


0. 


322 


1.614 




Avg. of Fits 


0.197 ±0.094 


0. 


323 ± 0.093 


1.610 ±0.088 




Theory 


0.200 ±0.082 


0. 


304 ±0.082 


1.535 ±0.082 


0.001-0.01 


Fit of Avg. 


0.202 


0. 


321 


1.605 




Avg. of Fits 


0.202 ±0.078 


0. 


321 ±0.077 


1.602 ±0.075 




Theory 


0.200 ±0.058 


0. 


303 ±0.058 


1.526 ±0.058 


0.002 -0.02 


Fit of Avg. 


0.206 


0. 


328 


1.594 




Avg. of Fits 


0.202 ±0.059 


0. 


324 ±0.057 


1.589 ±0.057 




Theory 


0.201 ± 0.041 


0. 


302 ± 0.041 


1.520 ±0.041 


0.0002 - 0.02 


Fit of Avg. 


0.206 


0. 


331 


1.609 




Avg. of Fits 


0.204 ±0.030 


0. 


330 ±0.031 


1.604 ±0.043 




Theory 


0.201 ±0.039 


0. 


303 ± 0.039 


1.527 ±0.039 



Table 2: PG-based fractal exponent estimates a for simulated FGNIF processes with a rate 
p = 40, for different frequency ranges. The governing rate processes have theoretical fractal 
exponents of a = 0.2, 0.8, and 1.5. The estimated fractal exponents were obtained from 
straight-line fits in doubly logarithmic coordinates to an average curve of 100 independent 
simulations (Fit of Avg.) and from averages of the slopes of the individual curves (Avg. 
of Fits). The deviations are the standard deviations obtained from the second averaging 
procedure. Also included are the theoretical bias values obtained from Eq. (32) and standard 
deviation values from [2, Eq. (17)]. These standard deviations are in substantial agreement 
with the data. Predictions of the bias underestimate the magnitude by factors of 3-10, 
although they correctly predict the sign. 

To demonstrate this, we estimated the fractal exponents over four time ranges (for the 
AF) and five frequency ranges (for the PG), with each range spanning a factor of ten or 
one hundred. As shown in Tables [l] and @, the results for the averaged AF and PG curves 
do not depend significantly on the choice of fit ranges, and in general agree well with the 
design values of the fractal exponents a. The deviations of the PG-based fractal exponent 
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estimator from these design values are accounted for in part by the theoretical bias values 
obtained from the general result provided in Eq. (|3~2"D , as shown in Table §; these theoretical 
predictions yield the correct sign for the bias, but underestimate the magnitude by factors 
of 3-10. Employing predictions directly based on the FGNIF process (see Appendix), rather 
than on a general FRSPP, yield results in substantial agreement with those of Eq. fl32|) . 
The AF-based fractal exponent estimators, in contrast, do not suffer from bias due to finite 
data size, and thus agree with the design values directly. In particular, the AF yields better 
performance than the PG estimate, especially for the larger design fractal exponents values 
of a = 0.8 and 1.5. 

The fractal exponents computed by the second method are also listed in Tables |1] and 0, 
indicated by "Avg. of Fit" , together with their standard deviations. (Note that for the first 
method there is no direct procedure to calculate the standard deviation since the results are 
extracted from a single curve). The two methods for estimating the fractal exponent agree 
remarkably well. For an ordinary Gaussian random variable, averaging over 100 samples 
would yield estimated values of the mean which exceeded the true mean by 2/10 of a standard 
deviation 5% of the time, with the rest clustered more tightly about the true mean. Of the 
27 experiments in which both methods for estimating the fractal exponent were performed, 
for only one (AF, a = 0.2, and 250 < T < 2500) does the sample mean exceed this limit, 
suggesting that this simple Gaussian model is valid. 

Useful estimators must have small variance as well as minimal bias; we therefore now focus 
on the standard deviations in Tables [I] and |2]. For both the AF- and PG-based estimators, and 
for ranges spanning one decade, the standard deviation decreases slightly as the design values 
of the fractal exponent increase; for two-decade ranges, this trend is reversed. Standard 
deviations for the two-decade ranges are significantly less than those for a single decade, as 
expected. For ranges spanning one decade, the standard deviation increase with decreasing 
frequency (PG) or increasing counting time (AF); at these time and frequency scales, the 
PG has the worst proportional frequency resolution, and the AF suffers from relatively few 
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counting windows over which to obtain accurate statistics. Finally, for similar time and 
frequency ranges the AF exhibits larger standard deviation than the PG. This is also likely 
due to the limited number of counting windows available, which proves more significant than 
the reduced resolution of the PG. Indeed, predictions of the estimated standard deviation for 
the PG based on earlier work || agree well with the data; these predictions in turn depend 
directly on the number of PG values available. 

A figure of merit for fractal exponent estimators which combines the bias and variance 
is the root-mean-square error, equivalent to the square root of the sum of the variance and 
the squared bias. Employing the results for the individual AF and PG entries provided in 
Tables |1] and |2] yields the results shown in Table |3|. For both the AF- and PG-based estimates, 



Statistic 


Range 


a = 0.2 


a = 0.8 


a = 1.5 


AF 


62.5 - 625 sec 


0.074 


0.072 


0.067 




125 - 1250 sec 


0.119 


0.114 


0.105 




250 - 2500 sec 


0.166 


0.153 


0.141 




25 - 2500 sec 


0.057 


0.058 


0.060 


PG 


0.00025 - 0.0025 Hz 


0.140 


0.142 


0.170 




0.0005- 0.005 Hz 


0.094 


0.096 


0.141 




0.001 - 0.01 Hz 


0.078 


0.080 


0.127 




0.002 - 0.02 Hz 


0.059 


0.062 


0.106 




0.0002 - 0.02 Hz 


0.030 


0.043 


0.113 



Table 3: Root-mean-square (rms) error of the estimated fractal exponent a for both AF- 
and PG-based estimators, for the same ranges and parameters used in Tables 1 and 2. 

the best performance is associated with times and frequencies well away from the duration 
of the data set: smallest lower counting time limits for the AF, and largest upper frequency 
limits for the PG, respectively. Comparing conjugate time and frequency ranges shows that 
the fractal-exponent estimator based on the AF yields significantly superior performance for 
the largest design value of a, while the PG proves somewhat better at the smallest value. 
Exact relationships between cutoff frequencies and cutoff times depend on the value of a 
itself [see also Eq. (|I8D1, so direct comparisons between ranges for the AF- and PG-based 
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estimates vary among the columns of Table [| 

The results obtained are governed by the finite duration of the data and by the nature 
of the estimation process, and are not an artifact of the particular value of the rate, p, that 
was chosen. Reducing the rate from p = 40 to p = 10, while keeping all other parameters 
fixed, generates results which differ from those presented in Tables |l| and @ by a maximum of 
±0.001. Thus in the FGNIF, the integrate-and-fire algorithm appears to translate the fractal 
quality of the continuous FGN into a point process nearly perfectly. All differences between 
the estimated fractal exponents and the design values used in the generation of the FGNIF 
simulations must therefore be due to finite-size effects: bias for the PG-based estimator, as 



discussed in Section B.2L and variance for both PG- and AF-based estimators. 



We now proceed to examine other point processes based on FGN, such as the FGNIF 
with jitter of the interevent intervals (FGNJIF) and the FGNDP. These processes, which 
include a second source of randomness that is inherent in the mechanism of point-process 
generation, exhibit increased fractal-exponent estimation error. We present averaged AF 
and PG plots for the FGNJIF with a mean rate p = 10 and jitter parameters a = 0, 0.5, 
and 1 in Fig. [8]. These plots attain asymptotic values at larger times (for the AF) and 
lower frequencies (PG) than those for the FGNIF do, especially for larger values of the jitter 
parameter a. This reduced range of scaling behavior serves to reduce the measured slopes 
of these curves and thus the estimated values of the fractal exponent for both the AF- and 
PG-based estimators, for all design values of a employed. For a — 0, corresponding to zero 
jitter (the FGNIF), the slopes attain the values shown in Tables [I] and [| Estimated fractal 
exponent values for the FGNJIF are presented in Table ^, which highlights the increasing 
negative bias for increased a. 

The AF curves in Fig. § exhibit two additional jitter-dependent features worthy of men- 
tion. First, as a increases the regularity of the signal decreases, leading to increased relative 
variance and a consequent decrease in the amplitude of the dip near T = 1 sec. Second, 
while the AF must approach an asymptotic value of unity for counting times T less than half 
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Statistic 


Range 


Jitter a 


a = 0.2 


a = 0.8 


a = 1.5 


AF 


250 - 2500 sec 





0.213 


0.804 


1.490 






0.5 


0.193 


0.796 


1.469 






1 


0.153 


0.779 


1.416 


PG 


0.002 - 0.02 Hz 





0.206 


0.828 


1.594 






0.5 


0.189 


0.796 


1.104 






1 


0.143 


0.718 


0.649 



Table 4: Fractal exponents estimated from the PG and AF of the FGNJIF process, with 
a chosen to be 0, 0.5, and 1. The governing rate processes have design fractal exponents of 
a = 0.2, 0.8, and 1.5, and mean values of p = 10 in all cases. One hundred individual curves 
were averaged, after which a single exponent was computed. As a increases the estimators 
exhibit an increasingly negative bias. 



the minimum interevent interval [see Eqs. flT3"| ) and (|T4])1, this does not appear to occur in 
Fig. H for <t = 1. This apparent conundrum is resolved by recognizing that large values of the 
jitter parameter a lead to an increased probability that two adjacent event times will be per- 
turbed to new values very close to each other. Thus as a increases, the relative frequency of 
quite small interevent intervals also increases, and the AF must be extended to ever smaller 
counting times to approach its limit of unity. We chose not to plot the AF curves at these 
smaller counting times T, since our focus lies in fractal behavior which manifests itself over 
longer time scales. 

We now turn to the FGNDP, which has a Poisson process substrate. As the average 
rate p decreases, the relative contribution of the randomness introduced by the Poisson 
generator increases as 1/^/p, and the PG- and AF-based fractal-exponents estimators exhibit 
increasingly larger bias. We present results for simulated FGNDP processes in Fig. ||; Table |5] 
provides estimated fractal exponent values. The average differences between the estimated 
values of the fractal-exponents a and the design values a are shown graphically in Fig. [10 
function of the mean event production rate p. 



42 



Statistic 


Range 


Rate p 


a = 0.2 


a = 0.8 


a = 1.5 


AF 


250 - 2500 sec 


5 


0.157 


0.774 


1.462 






20 


0.153 


0.794 


1.474 






40 


0.197 


0.798 


1.483 






80 


0.185 


0.802 


1.489 


PG 


0.002 - 0.02 Hz 


5 


0.130 


0.652 


0.997 






20 


0.132 


0.764 


1.092 






40 


0.149 


0.786 


1.351 






80 


0.157 


0.804 


1.496 



Table 5: Fractal exponents estimated from the PG and AF of the FGNDP process, with 
a rate p chosen to be 5, 20, 40, and 80. The governing rate processes have design fractal 
exponents of a = 0.2, 0.8, and 1.5. One hundred individual curves were averaged, after 
which a single exponent was computed. As p decreases the estimators exhibit an increasingly 
negative bias. 

For both the FGNJIF and the FGNDP processes, estimating fractal exponents over larger 
ranges of time or frequency leads to increased bias, particularly for large values of o and small 
values of p. Since the results depart from those obtained with the FGNIF, the additional 
randomness as p decreases, or as o increases, leads to a dilution of the power-law behavior of 
the AF and PG, changing its onset to at larger times and smaller frequencies. In the limits 
p — ► and a — > oo, the result is essentially white noise and the estimated fractal exponent 
will be zero. 

7 DISCUSSION 

Best results for the FGNIF were obtained with the AF-based fractal-exponent estimator 
using counting times in the range 25 < T < 2500; this yielded root-mean-square error 
results which were accurate to within 0.06 for all design values of a (Tables [l] and D). For 
both the AF- and PG-based estimators, performance degraded at longer times and smaller 
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frequencies (Tables [I], 0, and [3]). Finite-duration effects in the PG and relatively few counting 
windows in the AF produce significant error over these scales. The rms error remained near 
0.1 for both estimators, over all design values of the fractal exponent, and over all time and 
frequency ranges examined. 

The FGNJIF and FGNDP processes impose a second source of randomness on the fractal 
behavior of the FGN kernel. This additional randomness is essentially uncorrelated, and 
thus may be modeled as the linear superposition of white noise on the original fractal signal; 
this affects the shapes of both the AF and the PG. As this white-noise level increases, it 
dilutes the fractal structure present in the FGN process, beginning at the shortest time scales 
and moving progressively towards longer ones. This effect reduces the range of times and 
frequencies over which fractal behavior appears, resulting in a decreased estimated fractal 
exponent. This proved most significant for larger design values of a, and for the estimator 
based on the PG. Indeed, magnitudes of the bias exceeded 0.8 in one case. The AF-based 
estimator fared significantly better, although in this case it employed counting times which 
were five times the inverse of the frequencies used for the PG. These results compare favorably 
to earlier studies using the same number of events 0], for which inaccuracies of 0.1 or larger 
were obtained for a variety of fractal stochastic point processes, all of which included a 
second source of randomness. 

This study considered data sets with an average of 10 6 points; experimental data may 
contain far fewer points, and thus yield poorer fractal-exponent estimation statistics. In 
particular, finite-size effects, which lead to increased bias for the PG-based estimator and 
increased variance for the AF-based estimator will reduce the ranges of times and frequencies 
over which reliable estimation can be achieved. Similarly, data sets with larger numbers of 
points (computer network traffic, for example) would yield superior performance. In addition, 
if several independent data sets are available, known to be from identical experiments, then 
the AF plots or the resulting estimates could be averaged, reducing the effective variance. 
In this case, the reduced bias in the AF-based estimator would render it far superior to that 
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based on the PG. 

Finally we mention that there exist several methods for reducing the bias of the fractal 
exponent estimates if one has some a priori knowledge of the process of interest. Such an 
approach however does not follow the philosophy of estimating a completely unknown signal 
and we therefore did not attempt to compensate for the bias by the use of such methods. 

8 CONCLUSION 

We have investigated the properties of fractal and fractal-rate stochastic point processes 
(FSPPs and FRSPPs), focusing on the estimation of measures that reveal their fractal ex- 
ponents. The fractal-Gaussian-noise driven integrate-and-fire process (FGNIF) is unique as 
a point process in that it exhibits no short-term randomness; we have developed a gener- 
alization which includes jitter, the FGNJIF, for which the short-term randomness is fully 
adjustable. In addition to the randomness contributed by the fractal nature of the rate, all 
other FRSPPs appear to exhibit additional randomness associated with point generation, 
which confounds the accuracy of fractal-exponent estimation. The FGNJIF proved crucial 
in elucidating the role that such randomness plays in the estimation of fractal exponents 
in other FRSPPs. We presented analytical results for the expected biases of the PG- and 
AF-based fractional exponent estimators due to finite data length. We followed this theoret- 
ical work with a series of simulations of FRSPPs for three representative fractal exponents: 
a = 0.2, 0.8, and 1.5. Using these simulations together with the analytical predictions, 
we delineated the sources of error involved in estimating fractal exponents. In particular, 
using the FGNJIF we were able to separate those factors intrinsic to estimation over finite 
FRSPP data sets from those due to the particular form of FRSPP involved. We conclude 
that the AF-based estimate proves more reliable in estimating the fractal exponent than the 
PG-based method, yielding rms errors of 0.06 for data segments of FRSPPs with 10 6 points 
over all values of a examined. Finally, we note that wavelet generalizations of the AF appear 
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to yield comparable results |86|, ^7|], suggesting that the AF estimate may be optimal in some 
applications. 
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APPENDIX: PSD FOR THE FGNIF PROCESS 



Equation ( |32| ) yields predictions for the PG-based fractal exponent estimation bias for a 
general fractal-rate stochastic point process (FRSPP). More accurate results may be obtained 
for particular FRSPPs, taking into account the structure of these individual processes. We 
consider the case of the FGNIF process, with the further specialization that the original 
FGN array is obtained by direct Fourier synthesis with half of the array discarded, and that 
the analyzing periodogram subsequently doubles the number of points, yielding the same 
number as in the original array. 

We begin with a discrete-time power spectral density Sx(n) which is periodic with period 
M, and with Sx(0) = 0. For simulation purposes we define a conjugate symmetric sequence 
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X n based on the square root of Sx(n): 



X„ 







exp(jO n )JS x (n 



sgn(0 n 

Y"* 

^■M-n 



for n = 
for < n < M/2 
7r)y/S x (n) for n = M/2 

for M/2 <n< M, 



(35) 



where {#„} are independent and uniformly distributed in [0, 27r), sgn(-) represents the signum 
function, and * denotes complex conjugation. The corresponding time-domain signal may 
be obtained by inverse Fourier transformation: 
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We next define a subsampled sequence Z k , also of length M, by 
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for k even 
X( k -i)/2 for k odd, 



for which the corresponding sequence in the frequency domain becomes 
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Finally, we compute the expected value of the periodogram of Z 
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Numerical computation based on Eq. fl59| ) yields results substantially in agreement with 
those of Eq. ©. 
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Figure Captions 



Figure 1: Representations of a point process, (a) The events are represented by a sequence 
of idealized impulses, occurring at times t n , and forming a stochastic point process dN(t). 
For convenience of analysis, several alternative representations of the point process are used, 
(b) The counting process N(t). At every event occurrence the value of N(t) augments by 
unity, (c) The sequence of counts {Z k }, a discrete-time non- negative integer- valued stochas- 
tic process, is formed from the point process by recording the number of events in successive 
counting windows of length T. (d) The sequence of counts {Z k } can be conveniently de- 
scribed in terms of a count index k. Information is lost because the precise times of event 
occurrences within each counting window are eliminated in this representation. Correlations 
in the discrete-time sequence {Zk} can be readily interpreted in terms of real time. 

Figure 2: Doubly logarithmic plot of the count-based periodogram vs. frequency for the 
point process representing a spontaneous nerve spike train recorded at the output of a visual- 
system relay neuron in the lateral geniculate nucleus (LGN) of the cat (solid curve, cell No. 
MDS-LGN, X-ON type [63]). No external stimulus was present. The data segment analyzed 
here consists of 24285 events with an average interevent time of 0.133 sec, comprising a total 
duration of 3225 sec. Over long time scales (low frequencies), the curve can be fit by a straight 
line (dotted) of slope -0.67, representing fractal behavior. This least-squares straight-line fit 
to the data was calculated over the region from 0.002 Hz to 1 Hz. 

Figure 3: Doubly logarithmic plot of the Fano factor estimate vs. counting time, for the 
same spike train as used in Fig. 2 (solid curve). Over long time scales, the curve can be fit by 
a straight line (dotted) of slope 0.66, representing fractal behavior with a similar exponent 
to that obtained from the PG. This straight-line best fit to the data (on doubly logarithmic 
coordinates) was calculated over the region from 0.3 sec to 322 sec. 
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Figure 4: Doubly logarithmic plot of the Allan factor estimate vs. counting time for the 
same spike train as used in Figs. 2 and 3 (solid curve). Over long time scales, the curve can 
be fit by a straight line (dotted) of slope 0.65, representing fractal behavior with a similar 
exponent to that obtained from the PG and the FF. This straight-line best fit to the data 
(on doubly logarithmic coordinates) was calculated over the region from 0.3 sec to 322 sec. 

Figure 5: Sample functions of fractal renewal processes. Interevent times are power- law 
distributed, (a) The standard fractal renewal process (FRP) consists of Dirac 5 functions 
and is zero-valued elsewhere, (b) The alternating FRP switches between values of zero and 
unity. 

Figure 6: A sum of several independent and identical alternating FRPs (top) are added 
(center) to produce a fractal binomial noise process which serves as the rate function for 
an integrate-and-fire point process (bottom). The result is the fractal-binomial-noise-driven 
integrate-and-fire point process (FBNIF). 

Figure 7: A primary homogeneous Poisson point process M(t) with constant rate \i serves 
as the input to a linear filter with impulse response function h(t). The continuous-time 
stochastic process I(t) at the output of this filter is shot noise, which serves as the random 
rate for an integrate-and-fire process whose output is N(t). N(t) is the shot-noise driven 
integrate-and-fire point process (SNIF). If h(t) decays in a power-law fashion, then I(t) is 
fractal shot noise and N(t) is a fractal SNIF or FSNIF. 
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Figure 8: AF curves (left) and PG curves (right) of the fractal-Gaussian-noise driven 
integrate-and-fire process with jitter (FGNJIF), with a = 0, 0.5, and 1 for three design 
values of the fractal exponent: a = 0.2, 0.8, and 1.5. The mean rate p = 10 in all cases. 
The underlying rate function x[n] consisted of 2 15 samples, equally spaced at 1 sec each. 
For a = 0.2 and a = 0.8 the FF intercept time T was 1 sec; for a = 1.5 the PSD corner 
frequency u was 0.005 Hz. Results from 100 simulations were averaged. The jitter reduces 
the range of scaling and systematically lowers the average slopes of the curves. 



Figure 9: AF curves (left) and PG curves (right) of the fractal-Gaussian-noise driven 
Poisson process (FGNDP) with rate p = 5, 20, 40, and 80, for three design values of the 
fractal exponent: a = 0.2, 0.8, and 1.5. The underlying rate function x[n] consisted of 2 15 
samples, equally spaced at 1 sec each. For a = 0.2 and a = 0.8 the FF intercept time T 
was 10 times the average interevent time 1/p; for a — 1.5 the PSD corner frequency c^o was 
0.0005 times the average rate p. Results from 100 simulations were averaged. As the rate 
decreases, the bias of the PG- and AF-based fractal exponent estimators increases. 



Figure 10: Average differences between the estimated values of the fractal exponents a 
and the design values a = 0.2, 0.8, and 1.5 as a function of the mean event production rate 
p for the fractal-Gaussian- noise driven Poisson process (FGNDP). The AF-based estimate 
yields smaller bias than the PG-based measure. For both fractal exponent estimators, bias 
decreases with increasing average rate p, especially for the PG-based estimator with a = 1.5. 
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